*************************************************************
*** Monte Carlo experiments (SEM DGP) for the WWW JoP Project
***
*** Created: 9-19-14
*** Modified: 9-5-19
***
*************************************************************

clear
clear matrix
clear mata
set maxvar 32000

version 14

do "Programs\SAR Programs.do"

*************************************************************
*************************************************************
*** SEM Data-Generating Process
*************************************************************
*************************************************************


*************************************************************
*** Set up the parameters of the MC experiments
*************************************************************
scalar b1 = 1				/* Parameter for beta_x1 */
matrix I = I(200)			/* Identity matrix with N observations */

local run = 1				/* Run number (for running multiple do files at once) */
local s = 2					/* Standard deviation of the errors */
local MCsims = 1000			/* Number of simulations in Monte Carlo experiments */
local QIsims = 1000			/* Number of simulations for uncertainty in quantities of interest */
local obs = colsof(I)		/* Number of observations */
local seed = 648			/* Set the seed */

set seed `seed'
set obs `obs'
if `obs' > `QIsims' {
	set matsize `obs'
}
else {
	set matsize `QIsims'
}
*************************************************************

local st = "$S_TIME"
local oper = "SAR Model for `MCsims' simulations"

*** Load the n-order contiguity matrices
spatwmat using "Data\W\W_`obs'rs.dta", name(W)
spatwmat using "Data\W\W2_`obs'rs.dta", name(W2_c)
spatwmat using "Data\W\W3_`obs'rs.dta", name(W3_c)	

*** Generate second- and third-order W matrix
mat W2 = W*W
mat W3 = W*W*W

tempvar id
gen `id' = _n

preserve
	use "Data\W\W_`obs'rs.dta", clear
	gen id = _n
	spmat dta Wspat V*, id(id) replace
restore

mat cons = J(`obs',1,1)

*** Provide access to the matrices in mata
mata: I = st_matrix("I")
mata: W = st_matrix("W")
mata: W2 = st_matrix("W2")
mata: W3 = st_matrix("W3")
mata: W2_c = st_matrix("W2_c")
mata: W3_c = st_matrix("W3_c")
mata: beta = st_numscalar("b1")
local b1 = b1

local stfull = "$S_TIME"
local operfull = "SEM Model for `MCsims' simulations"

*************************************************************
*** Monte Carlo Experiments
*************************************************************
tempname sem
postfile `sem' sim lambda b1 str30 time b_x_sem se_x_sem lambda_sem se_lambda_sem teffect_sem deffect_sem ieffect_sem deffect_o0_sem deffect_o0_se_sem /* 
*/	true_teffect true_deffect true_ieffect true_feffect true_heffect true_deffect_o0 true_teffect_o1 true_deffect_o1 true_ieffect_o1 true_teffect_o2 true_deffect_o2 true_ieffect_o2 true_teffect_o3 true_deffect_o3 true_ieffect_o3 /*
*/	b_x_sar se_x_sar rho_sar se_rho_sar teffect_sar teffect_se_sar deffect_sar deffect_se_sar ieffect_sar ieffect_se_sar deffect_o0_sar deffect_o0_se_sar teffect_o1_sar teffect_o1_se_sar deffect_o1_sar deffect_o1_se_sar ieffect_o1_sar ieffect_o1_se_sar teffect_o2_sar teffect_o2_se_sar deffect_o2_sar deffect_o2_se_sar ieffect_o2_sar ieffect_o2_se_sar teffect_o3_sar teffect_o3_se_sar deffect_o3_sar deffect_o3_se_sar ieffect_o3_sar ieffect_o3_se_sar feffect_sar feffect_se_sar heffect_sar heffect_se_sar /*
*/	b_x_m1 se_x_m1 theta_1_m1 se_theta_1_m1 teffect_m1 teffect_se_m1 deffect_m1 deffect_se_m1 ieffect_m1 ieffect_se_m1 deffect_o0_m1 deffect_o0_se_m1 teffect_o1_m1 teffect_o1_se_m1 deffect_o1_m1 deffect_o1_se_m1 ieffect_o1_m1 ieffect_o1_se_m1 /*
*/	b_x_m2 se_x_m2 theta_1_m2 se_theta_1_m2 theta_2_m2 se_theta_2_m2 teffect_m2 teffect_se_m2 deffect_m2 deffect_se_m2 ieffect_m2 ieffect_se_m2 deffect_o0_m2 deffect_o0_se_m2 teffect_o1_m2 teffect_o1_se_m2 deffect_o1_m2 deffect_o1_se_m2 ieffect_o1_m2 ieffect_o1_se_m2 teffect_o2_m2 teffect_o2_se_m2 deffect_o2_m2 deffect_o2_se_m2 ieffect_o2_m2 ieffect_o2_se_m2 feffect_m2 feffect_se_m2 heffect_m2 heffect_se_m2  /*
*/	b_x_m3 se_x_m3 theta_1_m3 se_theta_1_m3 theta_2_m3 se_theta_2_m3 teffect_m3 teffect_se_m3 deffect_m3 deffect_se_m3 ieffect_m3 ieffect_se_m3 deffect_o0_m3 deffect_o0_se_m3 teffect_o1_m3 teffect_o1_se_m3 deffect_o1_m3 deffect_o1_se_m3 ieffect_o1_m3 ieffect_o1_se_m3 teffect_o2_m3 teffect_o2_se_m3 deffect_o2_m3 deffect_o2_se_m3 ieffect_o2_m3 ieffect_o2_se_m3 feffect_m3 feffect_se_m3 heffect_m3 heffect_se_m3  /*
*/	b_x_m4 se_x_m4 theta_1_m4 se_theta_1_m4 theta_2_m4 se_theta_2_m4 theta_3_m4 se_theta_3_m4 teffect_m4 teffect_se_m4 deffect_m4 deffect_se_m4 ieffect_m4 ieffect_se_m4 deffect_o0_m4 deffect_o0_se_m4 teffect_o1_m4 teffect_o1_se_m4 deffect_o1_m4 deffect_o1_se_m4 ieffect_o1_m4 ieffect_o1_se_m4 teffect_o2_m4 teffect_o2_se_m4 deffect_o2_m4 deffect_o2_se_m4 ieffect_o2_m4 ieffect_o2_se_m4 teffect_o3_m4 teffect_o3_se_m4 deffect_o3_m4 deffect_o3_se_m4 ieffect_o3_m4 ieffect_o3_se_m4 feffect_m4 feffect_se_m4 heffect_m4 heffect_se_m4  /*
*/	b_x_m5 se_x_m5 theta_1_m5 se_theta_1_m5 theta_2_m5 se_theta_2_m5 theta_3_m5 se_theta_3_m5 teffect_m5 teffect_se_m5 deffect_m5 deffect_se_m5 ieffect_m5 ieffect_se_m5 deffect_o0_m5 deffect_o0_se_m5 teffect_o1_m5 teffect_o1_se_m5 deffect_o1_m5 deffect_o1_se_m5 ieffect_o1_m5 ieffect_o1_se_m5 teffect_o2_m5 teffect_o2_se_m5 deffect_o2_m5 deffect_o2_se_m5 ieffect_o2_m5 ieffect_o2_se_m5 teffect_o3_m5 teffect_o3_se_m5 deffect_o3_m5 deffect_o3_se_m5 ieffect_o3_m5 ieffect_o3_se_m5 feffect_m5 feffect_se_m5 heffect_m5 heffect_se_m5  /*
*/	using "Generate\SEM\SEM.dta", replace every(20)

qui foreach r of numlist -0.8(.2)0.8 { 
	scalar lambda = `r'
	mata: lambda = st_numscalar("lambda")
	
	scalar true_teffect = `b1' + 0
	scalar true_deffect = `b1'
	scalar true_deffect_o0 = `b1'
	scalar true_ieffect = 0
	scalar true_feffect = 0
	scalar true_heffect = 0
	
	foreach o of numlist 1(1)3 {
		foreach e in t d i {
			scalar true_`e'effect_o`o' = 0
		}
	}		
	
	qui foreach i of numlist 1(1)`MCsims' {
		local st1 = "$S_TIME"
		local oper1 = "Simulation `i'"
		
		*** Generate the error term
		drawnorm e, n(`obs') sd(`s')
		mkmat e, matrix(e) 

		*** Generate the X matrix
		generate x1 = -10+(10 - -10)*uniform()
		mkmat x1, matrix(X)

		*** Generate the spatial lags
		mat slag1 = W * X
		svmat slag1 

		mat slag2 = W2 * X 
		svmat slag2

		mat slag3 = W3 * X
		svmat slag3

		mat slag2_c = W2_c * X 
		svmat slag2_c

		mat slag3_c = W3_c * X
		svmat slag3_c
		
		*=========================================================
		*=========================================================
		*** Generate the SEM data-generating process (mata)

		mata: X = st_matrix("X")
		mata: e = st_matrix("e")
		mata: y = (X*beta) + (lambda*W)*(e) + e // SEM model (Darmofal page 102; Anselin 2005 calls this a "spatial moving average")

		getmata y // extract the y from mata to view

		*============================================================
		*============================================================	
		if `obs' > `QIsims' {
			set obs `obs'
}
		else {
			set obs `QIsims'
		}		

		*** Estimate the SEM
		spreg ml y x1 in 1/`obs', id(`id') elmat(Wspat, eig)
		mat V = e(V)
		mat b = e(b)
		scalar b_x_sem = b[1,1]
		scalar se_x_sem = sqrt(V[1,1])
		
		scalar b_lambda_sem = _b[rho:_cons]			/* Rho and lambda are reversed in "spreg" Version 14 */
		scalar se_lambda_sem = _se[rho:_cons]
		
		scalar teffect_sem = b[1,1]
		scalar deffect_sem = b[1,1]
		scalar ieffect_sem = 0
		scalar deffect_o0_sem = b[1,1]
		scalar deffect_o0_se_sem = sqrt(V[1,1])
		
*		spreg_fit
*		scalar rmse_sem = r(rmse)


		
		*** Estimate the SAR
		spreg ml y x1 in 1/`obs', id(`id') dlmat(Wspat, eig)
		mat V = e(V)
		
		scalar b_x_sar = _b[y:x1]
		scalar se_x_sar = _se[y:x1]

		scalar b_rho_sar = _b[lambda:_cons]
		scalar se_rho_sar = _se[lambda:_cons]

		spreg_fit
		scalar rmse_sar = r(rmse)
		scalar rmse_xb_sar = r(rmse_xb)

		matrix coeff = [b_rho_sar,b_x_sar]
		matrix VCVM = [V[3,3], V[3,1] \ V[1,3], V[1,1]]
		
		preserve
			tempvar rho_draw beta_draw 
			corr2data `rho_draw' `beta_draw', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `rho_draw' `beta_draw', matrix(draws)
		restore		
		
		mata: draws = st_matrix("draws")
		mata: calceffsar(draws, I, W, `QIsims')
		
		scalar teffect_sar = r(teffect)
		scalar teffect_se_sar = r(teffect_se)
		scalar deffect_sar = r(deffect)
		scalar deffect_se_sar = r(deffect_se)
		scalar ieffect_sar = r(ieffect)
		scalar ieffect_se_sar = r(ieffect_se)
		
		scalar deffect_o0_sar = r(deffect0)
		scalar deffect_o0_se_sar = r(deffect0_se)
		
		scalar feffect_sar = r(feffect)
		scalar feffect_se_sar = r(feffect_se)
		scalar heffect_sar = r(heffect)
		scalar heffect_se_sar = r(heffect_se)
		
		foreach o of numlist 1(1)3 {
			foreach e in t d i {
				scalar `e'effect_o`o'_sar = r(`e'effect`o')
				scalar `e'effect_se_o`o'_sar = r(`e'effect`o'_se)
			}
		}
		
		* Model 1: W 
		qui reg y x1 slag11
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..2]
		matrix VCVM = V[1..2,1..2]
		
		preserve		
			tempvar beta_m1 theta1_m1
			corr2data `beta_m1' `theta1_m1', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m1' `theta1_m1', matrix(draws_m1)	
		restore		
		
		mata: draws_m1 = st_matrix("draws_m1")
		mata: calceffslx1(draws_m1, W, `QIsims')
	
		foreach e in t d i {	
			scalar `e'effect_m1 = r(`e'effect)
			scalar `e'effect_se_m1 = r(`e'effect_se)
			foreach o of numlist 1(1)1 {
				scalar `e'effect_o`o'_m1 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m1 = r(`e'effect_se_o`o')
			}			
		}
		
		scalar deffect_o0_m1 = r(deffect_o0)
		scalar deffect_se_o0_m1 = r(deffect_se_o0)
	
		scalar b_x_m1 = _b[x1]
		scalar se_x_m1 = _se[x1]
	
		scalar b_theta_1_m1 = _b[slag11]
		scalar se_theta_1_m1 = _se[slag11]

		
		
		* Model 2: W and W^2
		reg y x1 slag11 slag21
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..3]
		matrix VCVM = V[1..3,1..3]
		
		preserve		
			tempvar beta_m2 theta1_m2 theta2_m2
			corr2data `beta_m2' `theta1_m2' `theta2_m2', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m2' `theta1_m2' `theta2_m2', matrix(draws_m2)	
		restore		
		
		mata: draws_m2 = st_matrix("draws_m2")
		mata: calceffslx2(draws_m2, W, W2, `QIsims')

		foreach e in t d i {	
			scalar `e'effect_m2 = r(`e'effect)
			scalar `e'effect_se_m2 = r(`e'effect_se)
			foreach o of numlist 1(1)2 {
				scalar `e'effect_o`o'_m2 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m2 = r(`e'effect_se_o`o')
			}			
		}

		scalar feffect_m2 = r(feffect)
		scalar feffect_se_m2 = r(feffect_se)
		scalar heffect_m2 = r(heffect)
		scalar heffect_se_m2 = r(heffect_se)
		
		scalar deffect_o0_m2 = r(deffect_o0)
		scalar deffect_se_o0_m2 = r(deffect_se_o0)
		
		scalar b_x_m2 = _b[x1]
		scalar se_x_m2 = _se[x1]
	
		scalar b_theta_1_m2 = _b[slag11]
		scalar se_theta_1_m2 = _se[slag11]
	
		scalar b_theta_2_m2 = _b[slag21]
		scalar se_theta_2_m2 = _se[slag21]
	


		* Model 3: W and W2_c
		qui reg y x1 slag11 slag2_c1
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..3]
		matrix VCVM = V[1..3,1..3]
		
		preserve		
			tempvar beta_m3 theta1_m3 theta2_m3 
			corr2data `beta_m3' `theta1_m3' `theta2_m3', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m3' `theta1_m3' `theta2_m3', matrix(draws_m3)	
		restore		
		
		mata: draws_m3 = st_matrix("draws_m3")
		mata: calceffslx2(draws_m3, W, W2_c, `QIsims')

		foreach e in t d i {	
			scalar `e'effect_m3 = r(`e'effect)
			scalar `e'effect_se_m3 = r(`e'effect_se)
			foreach o of numlist 1(1)2 {
				scalar `e'effect_o`o'_m3 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m3 = r(`e'effect_se_o`o')
			}			
		}

		scalar feffect_m3 = r(feffect)
		scalar feffect_se_m3 = r(feffect_se)
		scalar heffect_m3 = r(heffect)
		scalar heffect_se_m3 = r(heffect_se)
		
		scalar deffect_o0_m3 = r(deffect_o0)
		scalar deffect_se_o0_m3 = r(deffect_se_o0)
		
		scalar b_x_m3 = _b[x1]
		scalar se_x_m3 = _se[x1]
	
		scalar b_theta_1_m3 = _b[slag11]
		scalar se_theta_1_m3 = _se[slag11]
	
		scalar b_theta_2_m3 = _b[slag2_c1]
		scalar se_theta_2_m3 = _se[slag2_c1]
	

		
		* Model 4: W, W^2 and W^3
		reg y x1 slag11 slag21 slag31
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..4]
		matrix VCVM = V[1..4,1..4]
		
		preserve		
			tempvar beta_m4 theta1_m4 theta2_m4 theta3_m4
			corr2data `beta_m4' `theta1_m4' `theta2_m4' `theta3_m4', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m4' `theta1_m4' `theta2_m4' `theta3_m4', matrix(draws_m4)	
		restore
		
		mata: draws_m4 = st_matrix("draws_m4")
		mata: calceffslx3(draws_m4, W, W2, W3, `QIsims')		
		
		foreach e in t d i {
			scalar `e'effect_m4 = r(`e'effect)
			scalar `e'effect_se_m4 = r(`e'effect_se)
						
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_m4 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m4 = r(`e'effect_se_o`o')
			}
		}		
		
		scalar feffect_m4 = r(feffect)
		scalar feffect_se_m4 = r(feffect_se)
		scalar heffect_m4 = r(heffect)
		scalar heffect_se_m4 = r(heffect_se)
		
		scalar deffect_o0_m4 = r(deffect_o0)
		scalar deffect_se_o0_m4 = r(deffect_se_o0)
		
		scalar b_x_m4 = _b[x1]
		scalar se_x_m4 = _se[x1]
	
		scalar b_theta_1_m4 = _b[slag11]
		scalar se_theta_1_m4 = _se[slag11]
	
		scalar b_theta_2_m4 = _b[slag21]
		scalar se_theta_2_m4 = _se[slag21]

		scalar b_theta_3_m4 = _b[slag31]
		scalar se_theta_3_m4 = _se[slag31]
		
		
		
		* Model 5: W, W2_c and W3_c
		qui reg y x1 slag11 slag2_c1 slag3_c1
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..4]
		matrix VCVM = V[1..4,1..4]
		
		preserve		
			tempvar beta_m5 theta1_m5 theta2_m5 theta3_m5
			corr2data `beta_m5' `theta1_m5' `theta2_m5' `theta3_m5', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m5' `theta1_m5' `theta2_m5' `theta3_m5', matrix(draws_m5)
		restore		
		
		mata: draws_m5 = st_matrix("draws_m5")
		mata: calceffslx3(draws_m5, W, W2_c, W3_c, `QIsims')

		foreach e in t d i {
			scalar `e'effect_m5 = r(`e'effect)
			scalar `e'effect_se_m5 = r(`e'effect_se)
						
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_m5 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m5 = r(`e'effect_se_o`o')
			}
		}		
		
		scalar feffect_m5 = r(feffect)
		scalar feffect_se_m5 = r(feffect_se)
		scalar heffect_m5 = r(heffect)
		scalar heffect_se_m5 = r(heffect_se)
		
		scalar deffect_o0_m5 = r(deffect_o0)
		scalar deffect_se_o0_m5 = r(deffect_se_o0)
		
		scalar b_x_m5 = _b[x1]
		scalar se_x_m5 = _se[x1]
	
		scalar b_theta_1_m5 = _b[slag11]
		scalar se_theta_1_m5 = _se[slag11]
	
		scalar b_theta_2_m5 = _b[slag2_c1]
		scalar se_theta_2_m5 = _se[slag2_c1]
	
		scalar b_theta_3_m5 = _b[slag3_c1]
		scalar se_theta_3_m5 = _se[slag3_c1]

		if mod(`i',5) == 0 {
			nois display "." _c
			if mod(`i',100) == 0 {
				nois display "$S_DATE : `c(current_time)': Simulation #`i'"
			}
		}		

		elapse "`st1'" "`oper1'"		
		post `sem' (`i') (`r') (`b1') ("$S_elap") (b_x_sem) (se_x_sem) (b_lambda_sem) (se_lambda_sem) (teffect_sem) (deffect_sem) (ieffect_sem) (deffect_o0_sem) (deffect_o0_se_sem) /* 		
		*/	(true_teffect) (true_deffect) (true_ieffect) (true_feffect) (true_heffect) (true_deffect_o0) (true_teffect_o1) (true_deffect_o1) (true_ieffect_o1) (true_teffect_o2) (true_deffect_o2) (true_ieffect_o2) (true_teffect_o3) (true_deffect_o3) (true_ieffect_o3) /*
		*/	(b_x_sar) (se_x_sar) (b_rho_sar) (se_rho_sar) (teffect_sar) (teffect_se_sar) (deffect_sar) (deffect_se_sar) (ieffect_sar) (ieffect_se_sar) (deffect_o0_sar) (deffect_o0_se_sar) (teffect_o1_sar) (teffect_se_o1_sar) (deffect_o1_sar) (deffect_se_o1_sar) (ieffect_o1_sar) (ieffect_se_o1_sar) (teffect_o2_sar) (teffect_se_o2_sar) (deffect_o2_sar) (deffect_se_o2_sar) (ieffect_o2_sar) (ieffect_se_o2_sar) (teffect_o3_sar) (teffect_se_o3_sar) (deffect_o3_sar) (deffect_se_o3_sar) (ieffect_o3_sar) (ieffect_se_o3_sar) (feffect_sar) (feffect_se_sar) (heffect_sar) (heffect_se_sar) /*
		*/	(b_x_m1) (se_x_m1) (b_theta_1_m1) (se_theta_1_m1) (teffect_m1) (teffect_se_m1) (deffect_m1) (deffect_se_m1) (ieffect_m1) (ieffect_se_m1) (deffect_o0_m1) (deffect_se_o0_m1) (teffect_o1_m1) (teffect_se_o1_m1) (deffect_o1_m1) (deffect_se_o1_m1) (ieffect_o1_m1) (ieffect_se_o1_m1) /*
		*/	(b_x_m2) (se_x_m2) (b_theta_1_m2) (se_theta_1_m2) (b_theta_2_m2) (se_theta_2_m2) (teffect_m2) (teffect_se_m2) (deffect_m2) (deffect_se_m2) (ieffect_m2) (ieffect_se_m2) (deffect_o0_m2) (deffect_se_o0_m2) (teffect_o1_m2) (teffect_se_o1_m2) (deffect_o1_m2) (deffect_se_o1_m2) (ieffect_o1_m2) (ieffect_se_o1_m2) (teffect_o2_m2) (teffect_se_o2_m2) (deffect_o2_m2) (deffect_se_o2_m2) (ieffect_o2_m2) (ieffect_se_o2_m2) (feffect_m2) (feffect_se_m2) (heffect_m2) (heffect_se_m2) /*
		*/	(b_x_m3) (se_x_m3) (b_theta_1_m3) (se_theta_1_m3) (b_theta_2_m3) (se_theta_2_m3) (teffect_m3) (teffect_se_m3) (deffect_m3) (deffect_se_m3) (ieffect_m3) (ieffect_se_m3) (deffect_o0_m3) (deffect_se_o0_m3) (teffect_o1_m3) (teffect_se_o1_m3) (deffect_o1_m3) (deffect_se_o1_m3) (ieffect_o1_m3) (ieffect_se_o1_m3) (teffect_o2_m3) (teffect_se_o2_m3) (deffect_o2_m3) (deffect_se_o2_m3) (ieffect_o2_m3) (ieffect_se_o2_m3) (feffect_m3) (feffect_se_m3) (heffect_m3) (heffect_se_m3) /*
		*/	(b_x_m4) (se_x_m4) (b_theta_1_m4) (se_theta_1_m4) (b_theta_2_m4) (se_theta_2_m4) (b_theta_3_m4) (se_theta_3_m4) (teffect_m4) (teffect_se_m4) (deffect_m4) (deffect_se_m4) (ieffect_m4) (ieffect_se_m4) (deffect_o0_m4) (deffect_se_o0_m4) (teffect_o1_m4) (teffect_se_o1_m4) (deffect_o1_m4) (deffect_se_o1_m4) (ieffect_o1_m4) (ieffect_se_o1_m4) (teffect_o2_m4) (teffect_se_o2_m4) (deffect_o2_m4) (deffect_se_o2_m4) (ieffect_o2_m4) (ieffect_se_o2_m4) (teffect_o3_m4) (teffect_se_o3_m4) (deffect_o3_m4) (deffect_se_o3_m4) (ieffect_o3_m4) (ieffect_se_o3_m4) (feffect_m4) (feffect_se_m4) (heffect_m4) (heffect_se_m4) /*
		*/	(b_x_m5) (se_x_m5) (b_theta_1_m5) (se_theta_1_m5) (b_theta_2_m5) (se_theta_2_m5) (b_theta_3_m5) (se_theta_3_m5) (teffect_m5) (teffect_se_m5) (deffect_m5) (deffect_se_m5) (ieffect_m5) (ieffect_se_m5) (deffect_o0_m5) (deffect_se_o0_m5) (teffect_o1_m5) (teffect_se_o1_m5) (deffect_o1_m5) (deffect_se_o1_m5) (ieffect_o1_m5) (ieffect_se_o1_m5) (teffect_o2_m5) (teffect_se_o2_m5) (deffect_o2_m5) (deffect_se_o2_m5) (ieffect_o2_m5) (ieffect_se_o2_m5) (teffect_o3_m5) (teffect_se_o3_m5) (deffect_o3_m5) (deffect_se_o3_m5) (ieffect_o3_m5) (ieffect_se_o3_m5) (feffect_m5) (feffect_se_m5) (heffect_m5) (heffect_se_m5)

		keep in 1/`obs'	
		drop e x1 y slag1* slag2* slag3*
		mata: mata drop draws draws_m1 draws_m2 draws_m3 draws_m4 draws_m5  /* fbeffect_m3 fbeffect_m5	*/
	}
}
	
postclose `sem'	
elapse "`stfull'" "`operfull'"
